program chem_gases_test

   use mod_chem
   use stdlib_string_type
   use iso_fortran_env, only: dp => real64
   implicit none

   type(type_chem_phases) :: phases
   integer, parameter :: n = 81
   real(dp) :: T(n), c(5, n), rho(n)
   integer :: file_io
   real(dp) :: x(n), y(n)
   integer :: j
   real(dp) :: tmp_i
   real(dp) :: rho_ref = 0.7915d-02, T_ref = 293.0_dp
   real(dp), allocatable :: production(:)

   phases = chem_file_load("air_5species_NASA9.yaml", "air")

   open(file_io, file = "gas_partial_i0.plt")
   do j = 1, n
      read(file_io, *)x(j), y(j), T(j), rho(j), c(1, j), c(2, j), c(3, j), c(4, j), c(5, j)
   enddo
   close(file_io)
   T = T*T_ref; rho = rho*rho_ref
   open(file_io, file = "gas_result.plt")

   write(file_io, *) 'VARIABLES = "X", "Y","wsi_O2", "wsi_N2",  "wsi_O",  "wsi_N",  "wsi_NO","mu", "k", "e"'
   ! T = 150.0_dp
   do j = 1, n
      call phases%set_state(T=T(j), Y=c(1:5, j), D=rho(j))
      production = phases%get_net_production_rates()*phases%get_molecular_weights()
      print*, j, T(j), rho(j), phases%get_internal_energy_mass()
      write(file_io, *) x(j)/1000.0_dp, y(j)/1000.0_dp, production(2), production(5), production(1), &
      production(3), production(4), phases%get_viscosity(T(j)), phases%get_thermal_conductivity(T(j)), &
      phases%get_internal_energy_mass(T(j))
   enddo
   close(file_io)

   open(file_io, file = "gamma.plt")
   write(file_io, *) 'VARIABLES = "X", "Y","gamma"'
   do j = 1, n
      call phases%set_state(T=T(j), Y=c(1:5, j), D=rho(j))
      print*, j, T(j), sqrt(phases%get_gamma()*phases%get_pressure()/phases%get_density()), &
      sqrt(phases%get_enthalpy_mass()/phases%get_internal_energy_mass()*phases%get_pressure()/phases%get_density())
      write(file_io, *) x(j)/1000.0_dp, y(j)/1000.0_dp, phases%get_gamma()
   end do
   close(file_io)

end
